1 Legende

Die folgende Notation wird in diesem Dokument benutzt:

ein_paar_R_befehle() #ein Kommentar dazu, wird von R ignoriert

## [1] "Die Ausgabe eines R-Befehls"

Aufgaben

Berlin

Das Arbeiten mit Geodaten ist eine zentrale Komponente der Umweltdatenverbeitung. Außerhalb Geographischer Informationssysteme (GIS) gibt es mächtige Werkzeuge, um direkt aus R heraus Geodaten zu analysieren und zu visualisieren. Im Gegensatz zur Erstellung einfacher Karten im GIS ist das zwar mit verhältnismäßig viel Code verbunden. Dafür hat man aber in R mehr Flexibilität und kann die Erstellung von Karten direkt in weitergehende Analysen einbetten. Auch die automatisierte Erstellung von Karten geht mit R leichter.

In diesem Kapitel werden wir Euch im Wesentlichen “Rezepte” vorstellen, wie Ihr unterschiedliche Typen von Geodaten in Eure R laden und damit arbeiten könnt. In diesem Kontext werden wir auch ein paar grundlegende Begriffe aus den Bereichen GIS und Kartographie rekapitulieren.

2 Karten erstellen in R…eigentlich gar nichts Neues?

Ihr habt in diesem Kurs bereits einge Pakete und Funktionen für die visuelle Darstellung von Daten kennengerlernt. Diese Pakete ermöglichen es, Punkte, Linien oder Flächen in ein Koordinationsystem zu plotten - und um nichts anderes handelt es sich bei einer Karte.

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

Für das folgende Beispiel benötigen wir zunächst das Paket sf. Installiere es, falls nötig.

library(sf)

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

Wir lesen die Metadaten der Niederschlagsstationen des DWD ein. Nettes feature: Wir müssen die Daten gar nicht runterladen, sondern können direkt die Datei via https ansprechen.

fdwdgauges <- "https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/daily/more_precip/recent/RR_Tageswerte_Beschreibung_Stationen.txt"
dwdgauges <- read.fwf(fdwdgauges, c(5,9,9,15,12,10, 42, 201), 
                     encoding="ISO-8859-1", skip=2,header = FALSE,
                     col.names =  c("id", "from", "to", "rioation", 
                                    "lat", "lon","name","state"))
head(dwdgauges)

##   id     from       to rioation     lat     lon
## 1  1 19120101 19860630      478 47.8413  8.8493
## 2  2 19510101 20061231      138 50.8066  6.0996
## 3  3 18910101 20110331      202 50.7827  6.0941
## 4  4 19510101 19791031      243 50.7683  6.1207
## 5  6 19821101 20210927      455 48.8361 10.0598
## 6  7 19510101 19960131      473 48.8140 10.1285
##                                         name
## 1  Aach                                     
## 2  Aachen (Kläranlage)                      
## 3  Aachen                                   
## 4  Aachen-Brand                             
## 5  Aalen-Unterrombach                       
## 6  Aalen-Unterkochen                        
##                                                                                                state
## 1 Baden-Württemberg                                                                                 
## 2 Nordrhein-Westfalen                                                                               
## 3 Nordrhein-Westfalen                                                                               
## 4 Nordrhein-Westfalen                                                                               
## 5 Baden-Württemberg                                                                                 
## 6 Baden-Württemberg

Wir sehen: Die Datei enthält die geographischen Koordinaten der Stationen (lon, lat). Damit können wir die Positionen auf eine Karte plotten.

plot(lat ~ lon, dwdgauges, asp = 1,
     main = "DWD Niederschlagsstationen")

Wow…ganz schön viele Stationen! Und die messen alle Niederschlag?!

Manipulieren wir doch mal die Karte, indem wir z.B. diejenigen Stationen in rot darstellen, die gar nicht mehr aktiv sind. Dazu müssen wir die Zeitkomponente der Metadaten als Datum interpretieren und anschließend den dataframe mit Hilfe des to-Attributs filtern.

dwdgauges$fromdate <- as.Date(as.character(dwdgauges$from), format="%Y%m%d")
dwdgauges$todate <- as.Date(as.character(dwdgauges$to), format="%Y%m%d")
inactive <- dwdgauges$todate < "2021-07-01"

plot(lat ~ lon, dwdgauges[!inactive,], asp = 1,
     main = "DWD Niederschlagsstationen")
points(lat ~ lon, dwdgauges[inactive,], asp = 1, col = "red")

Aha…also doch nicht so viele? Aber klar: Stationen wurden über die Zeit verlegt, stillgelegt, neu gebaut und wieder stillgelegt. Wenn man die Daten verarbeiten möchte, muss aber jeder Standort eine eindeutige ID haben…da kommt dann eben über die Zeit was zusammen. Apropos Zeit…reichern wir doch mal die Karte mit einer Zeitreihe an: Zahl der Stationen von 1901 bis 2021.

Installiere und lade das Paket lubridate, bevor Du die nächsten Zeilen ausführst. Das macht das Arbeiten mit Datumsdaten leichter.

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

# Heute ist welches Jahr?
endyear   <- year(today())

# Vektor mit Jahreszahlen von 1901-heute
years     <- 1901:endyear

# Preallokation eines Vektors in den wir das Zählergebnis für jedes 
# Jahr reinschreiben
numgauges <- rep(0,length(years))

# Um uns die Arbeit zu erleichtern 
# Man erinnere sich: Stay DRY, don’t get WET!
hasdatafun <- function(data,Y){
  hasdata = Y >= year(data$fromdate) & 
            Y <= year(data$todate)
}

# Anhand einer Schleife zählen wir nun für jedes Jahr die Anzahl der 
# Stationen
for (iter in 1:length(years)){
  hasdata = hasdatafun(dwdgauges,years[iter])
  numgauges[iter] <- sum(hasdata)
}

# Nun können wir die Daten plotten
par(mfrow = c(1,4))
hasdata <- hasdatafun(dwdgauges,1920)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = "1920", pch = 20)

hasdata <- hasdatafun(dwdgauges,1987)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = "1987",ylab = "", pch = 20)

hasdata <- hasdatafun(dwdgauges,endyear)
plot(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = as.character(endyear),ylab = "", pch = 20)

plot(years,numgauges,type = "l")

Fazit

Wir haben einen Geodatensatz eingelesen und Karten erstellt, ohne ein neues Werkzeug zu nutzen. Obendrein haben wir die Karten ganz lässig mit Zeitreihen angereichert… an sowas müsstet Ihr in einem GIS lange schrauben!

Trotzdem sind wir noch nicht zufrieden. Es wäre z.B. schön, die Landesgrenzen und ggf. Grenzen der Bundesländer mit darzustellen…

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

Wir haben nun also mit Standardpaketen bereits wesentliche Elemente von Karten erstellt? Was fehlt uns noch?

Einlesen gängiger Geodaten in Vektor- und Rasterformat.
Koordinatentransformation
Effizientes Plotten spezifischer Geodatentypen- und Elemente.
Gängige GIS-Operationen

4 Aufgabe

Diskutiert (am Besten mit Euren Nachbar:innen) mal ein paar GIS-Grundbegriffe…

Was sind Vektordaten? Nennt zwei gängige Datenformate!
Was sind Rasterdaten? Nennt zwei gängige Datenformate!
Nennt mindestens drei typische GIS-Operationen mit Vektordaten.
Wozu braucht man Projektionen (oder Spatial Reference Systems)?

5 Lesen und plotten von Vektordaten mit `sf`

Die folgende Abbildung macht noch einmal deutlich, worum es sich bei Vektordaten handelt.

Vektordaten

(Quelle: Colin Williams, NEON)

Typische Datenformate für Vektordaten:

Shapefile (.shp),
geojson (.json)
GPS Exchange Format (.gpx)
Google Keyhole Markup Language (.KML, .KMZ)
…

6 Aufgabe

Diskutiert mit Euren Nachbarn: Für welche Umweltdaten ist das Vektorformat geeeinget, für welche nicht?

Vektordaten sind geeignet für:

Verwaltungeinheiten, politische Grenzen

Straßen, andere Infrastrukturelemente

Gebäudeumrisse

Gewässer (Flüsse, Seen, …)

Habitate

Messpunkte

…

Vektordaten sind ungeignet für:

Digitale Geländemodelle

Fernerkundungsdaten

Daten aus Klimamodellen

…

7 Vektordaten mit `sf`

sf hält wesentliche Funnktionen bereit, die wir uns in Python für Vektordaten wünschen.

Beginnen wir doch damit, unserer Karte Landesgrenzen hinzuzufügen. Das Shapefile mit den Bundesländergrenzen haben wir direkt vom Bundesamt für Kartographie und Geodäsie heruntergeladen (Link).

library(sf)

f_bundeslaender <- "../data/shapefiles/VG250_LAN.shp"
bundeslaender   <- st_read(f_bundeslaender)

## Reading layer `VG250_LAN' from data source 
##   `C:\Users\Wolf\ownCloud\Lehre\GEE-M-MK6_Datenverarbeitung\umweltdv\06_geodata_maps\data\shapefiles\VG250_LAN.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 35 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 280371.1 ymin: 5235856 xmax: 921292.4 ymax: 6106244
## Projected CRS: ETRS89 / UTM zone 32N

head(bundeslaender)

## Simple feature collection with 6 features and 26 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 280371.1 ymin: 5471422 xmax: 674202.3 ymax: 6101444
## Projected CRS: ETRS89 / UTM zone 32N
##   ADE GF BSG ARS AGS      SDV_ARS                 GEN                  BEZ IBZ
## 1   2  4   1  01  01 010020000000  Schleswig-Holstein                 Land  20
## 2   2  4   1  02  02 020000000000             Hamburg Freie und Hansestadt  22
## 3   2  4   1  03  03 032410001001       Niedersachsen                 Land  20
## 4   2  4   1  04  04 040110000000              Bremen     Freie Hansestadt  23
## 5   2  4   1  05  05 051110000000 Nordrhein-Westfalen                 Land  20
## 6   2  4   1  06  06 064140000000              Hessen                 Land  20
##   BEM NBD SN_L SN_R SN_K SN_V1 SN_V2 SN_G FK_S3 NUTS        ARS_0    AGS_0
## 1  --  ja   01    0   00    00    00  000     0  DEF 010000000000 01000000
## 2  --  ja   02    0   00    00    00  000     0  DE6 020000000000 02000000
## 3  --  ja   03    0   00    00    00  000     0  DE9 030000000000 03000000
## 4  --  ja   04    0   00    00    00  000     0  DE5 040000000000 04000000
## 5  --  ja   05    0   00    00    00  000     0  DEA 050000000000 05000000
## 6  --  ja   06    0   00    00    00  000     0  DE7 060000000000 06000000
##          WSK         DEBKG_ID RS       SDV_RS         RS_0
## 1 2012-02-01 DEBKGDL200000029 01 010020000000 010000000000
## 2 1974-01-01 DEBKGDL20000E6GD 02 020000000000 020000000000
## 3 2015-01-01 DEBKGDL20000E6EW 03 032410001001 030000000000
## 4 2010-01-01 DEBKGDL20000E0SF 04 040110000000 040000000000
## 5 2009-11-01 DEBKGDL20000E6GR 05 051110000000 050000000000
## 6 2015-01-01 DEBKGDL20000E3LK 06 064140000000 060000000000
##                         geometry
## 1 MULTIPOLYGON (((464810.7 61...
## 2 MULTIPOLYGON (((578219 5954...
## 3 MULTIPOLYGON (((479451.1 59...
## 4 MULTIPOLYGON (((466930.3 58...
## 5 MULTIPOLYGON (((477607.3 58...
## 6 MULTIPOLYGON (((534242 5721...

Puuh, 26 Spalten - ganz schön viele Attribute…naja, wir intessieren uns ja erstmal nur für die Ländergrenzen. Schnell mal gucken, ob die da sind…

Die Funktion st_geometry legt erstmal alle Attribute beiseite, so dass nur die Geometrie geplottet wird.

plot(st_geometry(bundeslaender), axes = T)

Ok, sieht ja schon gut aus…sogar mit Hoheitsgebieten im Meer - wollten wir die haben? Egal.

Was aber bei den Achsen auffällt, ist die Größenordnung…irgendwas um 600.000 für die x-Achse und zwischen 5 und 6 Millionen für die y-Achse. Längen- und Breitengrade sind das jedenfalls nicht.

st_crs(bundeslaender)

## Coordinate Reference System:
##   User input: ETRS89 / UTM zone 32N 
##   wkt:
## PROJCRS["ETRS89 / UTM zone 32N",
##     BASEGEOGCRS["ETRS89",
##         DATUM["European Terrestrial Reference System 1989",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4258]],
##     CONVERSION["UTM zone 32N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Europe between 6Â°E and 12Â°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
##         BBOX[38.76,6,83.92,12]],
##     ID["EPSG",25832]]

Bevor wir Euch ins offene Messer laufen lassen: Wir müssen also umprojizieren. Die Daten liegen im projezierten Koordinatensystem UTM Zone 32 vor. Für unseren Zweck ist ein geographisches Koordinatensystem besser. Die EPSG Nummer 4326 sollte man sich merken. Sie steht für das WGS84 Koordinatensystem.

bundeslaender2 <- st_transform(bundeslaender, crs = 4326)

Nun können wir unsere Stationsdaten zusammen mit den Bundesländern plotten.

plot(st_geometry(bundeslaender2), axes = T, main = "Stationen im Jahr 1920")

hasdata <- hasdatafun(dwdgauges,1920)
points(lat ~ lon, dwdgauges[hasdata,], asp = 1,
     main = as.character(endyear))

Oder wir projizieren die Stationsdaten um. Dazu müssen wir den Dataframe zunächst in ein sf-Objekt umwandeln.

dwdgaugessf <- st_as_sf(dwdgauges, coords = c("lon","lat"), crs = 4326)
dwdgaugessf <- st_transform(dwdgaugessf, crs = st_crs(bundeslaender))

plot(st_geometry(bundeslaender), axes = T, main = "Stationen im Jahr 1920")

hasdata <- hasdatafun(dwdgauges,1920)
plot(st_geometry(dwdgaugessf[hasdata,]),
     main = as.character(endyear), add = T)

8 Rasterdaten mit `raster`

raster ist ein Paket zum Lesen, Schreiben und Verarbeiten von Rasterdaten. Eine ausführliche Dokumentation des Pakets findet Ihr hier. Ein weiteres, häufig verwendetes Paket sind stars (hier) und terra (hier).

Wir werden raster in diesem Kurs v.a. zum Lesen eines Rasters verwenden. Unser Ziel ist es, unsere schöne Deutschlandkarte mit einer Topographie zu hinterlegen. Es gibt für die meisten Bundesländer mittlerweile frei verfügbare hochaufgelöste (1x1 m) Digitale Geländemodelle (DGM; engl. digital rioation model, DEM) . Außerdem gibt es hochauflösende globale DGM mit Auslösungen zwischen 20 und 30 Metern (u.a. SRTM, ASTER). Eine “optimale” Kombination von SRTM und ASTER auf europäischer Ebene gibt es bei der Europäische Umweltagentur EEA hier.

Für eine deutschlandweite Darstellung sind diese Produkte aber zu hoch aufgelöst und damit ineffizient (Speicher, Rechenbedarf). Wir nutzen daher ein Produkt, das auf eine Auflösung von 1x1 km für ganz Europa “resampled” wurde. Es beruht auf dem globalen Datensatz GTOPO30. Auch diesen Datensatz haben wir bei der EEA hier heruntergeladen.

library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/Wolf/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Wolf/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/Wolf/Documents/R/win-library/4.1/rgdal/proj

library(raster)

ELEV <- raster('../data/raster/elevation1x1_new.tif')

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps User_Defined_Spheroid in Proj4 definition:
## +proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +R=6378137 +units=m
## +no_defs +type=crs

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum User_Defined in Proj4 definition

Das erhaltene Objekt rio ist vom Typ raster und wir können es ein wenig interviewen, ohne überhaupt etwas in den Arbeitsspeicher geladen zu haben.

cat(paste("Beginn des Interviews am", now(),'\n'))

## Beginn des Interviews am 2021-09-29 09:01:58

cat(paste("Wie heißen Sie?:\n", filename(ELEV),'\n'))

## Wie heißen Sie?:
##  C:\Users\Wolf\ownCloud\Lehre\GEE-M-MK6_Datenverarbeitung\umweltdv\06_geodata_maps\data\raster\elevation1x1_new.tif

cat(paste("Wie viele Zeilen haben Sie?\n", nrow(ELEV),'\n'))

## Wie viele Zeilen haben Sie?
##  4251

cat(paste("Wie viele Spalten haben Sie?\n", ncol(ELEV),'\n'))

## Wie viele Spalten haben Sie?
##  4901

cat(paste("Wie lautet Ihre Bounding Box?\n", extent(ELEV),'\n'))

## Wie lautet Ihre Bounding Box?
##  extent(2399843.1035138, 7300843.1035138, 1249900.3833563, 5500900.3833563)

Vor allem aber die Frage, die wir jedem Geodatensatz stellen müssen: Welche Projektion (“coordinate reference system”, CRS)?

cat(paste("Welches Koordinatensystem haben Sie?\n", crs(ELEV),'\n'))

## Welches Koordinatensystem haben Sie?
##  +proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +R=6378137 +units=m +no_defs

Wir haben es also mit laea, der Lambert Azimuthal Equal Area, also einer flächentreuen Projektion, zu tun. Nun gut, wir nehmen das einfach mal so hin.

Rasterdatenformate wie GeoTIFF sind in der Fernerkundung beliebt. Daher sind sie so ausgelegt, dass in einem Datensatz unterschiedliche Bänder (engl. “bands”) abgelegt werden können (klassischerweise z.B. rot, grün, blau - RGB). Im unserem Datensatz erwarten wir aber nur die Geländehöhe. Dies können wir verifizieren:

cat(paste("Wie viele Bänder haben Sie?\n", nbands(ELEV),'\n'))

## Wie viele Bänder haben Sie?
##  1

Ok, werfen wir mal einen Blick auf die Werte, die in dem Raster gespeichert sind. Fragen wir doch nach dem größten Wert.

cat(paste("Welches ist ihr maximaler Wert?\n", maxValue(ELEV),'\n'))

## Welches ist ihr maximaler Wert?
##  187

9 Aufgabe

Stopp…die maximale Geländehöhe in Europa soll 187 Meter betragen? Überlegen Sie, was für ein Missverständnis/Problem hier vorliegen könnte. Lesen Sie ggf. die Metadaten (Link). Schlagen Sie eine pragmatische Lösung vor.

In den Metadaten finden wir, dass ein z factor für den Datensatz genutzt wurde. Wir müssen also das Raster mit dem Wert 10 multiplizieren.
ELEV <- ELEV*10

10 Testplot

plot(ELEV)

Soweit, so gut: Wir haben einen Rasterdatensatz mit raster eingelesen, wir kennen sein CRS und können den Datensatz mit plot plotten.

11 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

Wir plotten alles in das CRS des Rasterdatensatzes, weil wir uns mit dem Thema Projektion bzw. Warpen von Rasterdaten nicht befassen wollen. Wir bringen also zunächst unsere Vektordaten (Bundesländer und Stationsmessnetz) in das CRS des Rasterdatensatzes.

bundeslaender_laea <- st_transform(bundeslaender, crs = crs(ELEV))
dwdgaugessf  <- st_transform(dwdgaugessf, crs = crs(ELEV))

Und weil die Topographie mit einem Hillshade immer schicker aussieht, nutzen wir noch das Hillshade-Raster, welches uns auf der gleich EEA-Seite zur Verfügung gestellt wird (Link). Das Hillshade-Raster plotten wir ebenfalls mit imshow, aber mit einer Tansparenz (Parameter alpha). Das Ganze machen wir zudem mit dem Paket tmap, einem Paket für die Erstellung thematischer Karten.

hillsh = raster('../data/raster/hillshade1x1/hillshade1x1.tif')

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps User_Defined_Spheroid in Proj4 definition:
## +proj=laea +lat_0=52 +lon_0=20 +x_0=5071000 +y_0=3210000 +R=6378137 +units=m
## +no_defs +type=crs

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum D_User_Defined in Proj4 definition

library(tmap)

Da wir nicht ganz Europa plotten wollen, geht es schneller, wenn wir die beiden Rasterdaten “clippen”.

hillshc <- crop(hillsh,extent(bundeslaender_laea))/255
ELEVc <- crop(ELEV,extent(bundeslaender_laea))
ELEVc[ELEVc == 0] = NA

tmap_mode("plot")

## tmap mode set to plotting

## tmap mode set to plotting
tm_shape(ELEVc) +
    tm_raster("elevation1x1_new", palette = terrain.colors(255), 
              legend.show = FALSE) +
tm_shape(hillshc) +
    tm_raster("hillshade1x1", palette = gray.colors(255), 
              alpha = 0.5, legend.show = FALSE) +  
tm_shape(bundeslaender_laea) +
    tm_borders("white", lwd = .5) +
tm_shape(dwdgaugessf[hasdata,]) +
    tm_symbols(col = "black", size = 0.3) +
tm_legend(show = FALSE) +
tm_compass() + 
tm_graticules() +
tm_layout(title = "DWD Niederschlagsstationen 1920", 
          title.bg.color = "white", title.bg.alpha = 0.7)

Alternativ können wir mit tmap auch eine interaktive Webkarte erstellen.

tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(dwdgaugessf[hasdata,]) +
    tm_symbols(col = "black", size = 0.3)

12 Übung für die Coding-Werkstatt

Fertige auf Basis der heute behandelten Inhalte und Daten eine Karte für Dein Heimatbundesland an, in welchem folgende Informationslayer vorhanden sind:

Topographie
Landesgrenzen
Hauptstadt
wesentliche Fließgewässer
…

Sucht Euch ggf. erforderliche Geodaten selbst aus dem Netz!

13 Weiteres Studium

Viele weitere Infos und Beispiele findet ihr auf der Website RSpatial und r-spatial.

Fine

Projektionen

Quelle: https://xkcd.com/977/

Geodaten in R - mehr als nur eine Alternative zum GIS

Eine Karte sagt mehr als tausend Worte

Wolfgang Schwanghart

1 Legende

Aufgaben

2 Karten erstellen in R…eigentlich gar nichts Neues?

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

4 Aufgabe

5 Lesen und plotten von Vektordaten mit `sf`

6 Aufgabe

7 Vektordaten mit `sf`

8 Rasterdaten mit `raster`

9 Aufgabe

10 Testplot

11 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

12 Übung für die Coding-Werkstatt

13 Weiteres Studium

Geodaten in R - mehr als nur eine Alternative zum GIS

Eine Karte sagt mehr als tausend Worte

Wolfgang Schwanghart

1 Legende

Aufgaben

2 Karten erstellen in R…eigentlich gar nichts Neues?

2.1 Beispiel: Niederschlagsstationskollektiv des Deutschen Wetterdienstes

3 Was fehlt uns, um effizient schöne Karten zu gestalten?

4 Aufgabe

5 Lesen und plotten von Vektordaten mit sf

6 Aufgabe

7 Vektordaten mit sf

8 Rasterdaten mit raster

9 Aufgabe

10 Testplot

11 Synthese: Stationsmessnetz mit Ländergrenzen und Topographie plotten

12 Übung für die Coding-Werkstatt

13 Weiteres Studium

5 Lesen und plotten von Vektordaten mit `sf`

7 Vektordaten mit `sf`

8 Rasterdaten mit `raster`